Introducción a R

Probabilidad y Estadística Fundamental

Lenguaje R

R provee un rico ambiente para el análisis de datos y gráficos. Este puede ser descargado libremente desde el sitio https://cran.r-project.org/ para Windows, Linux y Mac. La instalación es sencilla y en este momento se asume que el programa ya ha sido descargado e instalado en su computador.

Objetivo: Introducir lenguaje R enfoncandose principalmente en los items que son útiles para el desarrollo de actividades de la asignatura Probabilidad y Estadística Fundamental. Tales como

• Léxico básico (variables, constantes, paquetes, directorio de trabajo)
• Operaciones básicas
• Operaciones aritméticos
• Estructuras de datos
• Gráficos y medidas descriptivas univariadas
• Gráficos y medidas descriptivas multivariadas
• Representación tabular de variables
• Funciones de probabilidad y de densidad (discretas y continuas)
• Calculo de cuantiles de la distibución normal y t-student.
• Estimación por intervalos y pruebas de hipótesis.
• Estimación de parámetros de un modelo de regresión.

R Studio

La programación se realizará sobre ambiente de desarrollo RStudio que se puede descargar de manera libre en https://www.rstudio.com/products/rstudio/download/ previa instalación de R.

Léxico básico

Expresiones aritméticas

Después de iniciar R, iniciemos con algunos comandos en la consola, por ejemplo, con algunas expresiones aritméticas simples

2+1   # Suma

[1] 3

2*3   # Multiplicación

[1] 6

2^2   # Potencia

[1] 4

2**3  # Potencia

[1] 8

6/2   # División

[1] 3

Note que es posible incluir comentarios en el código usando el simbolo #

Funciones aritméticas

También todas las funciones aritméticas tal como \(\log(), \exp(), \sin(),\pi ,\sqrt()\) están disponibles

log(exp(1))

[1] 1

log(10, base = 10)

[1] 1

exp(1)

[1] 2.718282

sin(0)

[1] 0

pi

[1] 3.141593

sqrt(25)

[1] 5

Asignar un valor a una variable

Para la asignación de un valor a una variable (que usualmente se representa mediante letras) se usará cualquiera de los operadores =, <-. Se imprimirá el valor de la variable mediante la función print() aunque llamando la variable también retornará el valor.

x = 2
print(x)

[1] 2

y <- 3+5; print(y)

[1] 8

print(x+y^2)

[1] 66

Nota: se podrán escribir en la misma línea del código varias intrucciones, recuerde separar cada una mediante punto y coma ;.

Operadores lógicos en R

Objetos de R

Vectores

Para crear vectores es usada la función c() (concatenar o coleccionar), se le asigna el vector a una variable del mismo modo.

a = c(2,4,6,8,10); print(a)

[1]  2  4  6  8 10

s = c("A","B","C","A","C","D","A"); print(s)

[1] "A" "B" "C" "A" "C" "D" "A"

Nota: Se pueden crear vectores de caracteres tal como el vector s.

Se puede comprobar el tipo de objeto de a y s con la función class().

class(a)   # vector numérico

[1] "numeric"

class(s)   # vector de caracteres

[1] "character"

Operaciones con vectores

En primer lugar, obtenemos la longitud de los vectores mediante la función length()

length(a)

[1] 5

length(s)

[1] 7

Operaciones aritméticas pueden llevarse a cabo, así si se quiere multiplicar el vector por una constante (escalar) o por otro vector el resultado será la multiplicaciónd ela constante por cada componente del vector.

a2 = a*3
print(a); print(a2)

[1]  2  4  6  8 10

[1]  6 12 18 24 30

Del mismo modo si se desea utilizar funciones aritméticas. Estas afectarán cada componente del vector.

a3 = sqrt(a)
print(a); print(a3)

[1]  2  4  6  8 10

[1] 1.414214 2.000000 2.449490 2.828427 3.162278

Las operaciones también se podrán realizar entre dos o más vectores.

a4 = a + a2^2 + log(a3)
print(a4)

[1]  38.34657 148.69315 330.89588 585.03972 911.15129

Importante: Las operaciones entre vectores se hacen con aquellos de misma longitud. Note que en los ejemplos realizados se están operando vectores de longitud 5.

Vector Secuencia o repetición

En muchos casos requerimos vectores que contengan una progresión numérica, o un vector cuyas componentes son las mismas, para esto R dispone de las funciones seq() y rep(). Observe los siguientes ejemplos.

s1 = seq(from = 2, to = 27) ## vector definido de 2 a 27
print(s1); length(s1)

 [1]  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
[24] 25 26 27

[1] 26

s2 = seq(from =2, to = 27, length.out = 6)  ## Longitud de salida de 6, secuencia cada 5
print(s2); length(s2)

[1]  2  7 12 17 22 27

[1] 6

s3 = seq(from =2, to = 27, by = 3)  ## Longitud de salida de 9, secuencia cada 3
print(s3); length(s3)

[1]  2  5  8 11 14 17 20 23 26

[1] 9

Nota: from=, to=, lenght.out=, by= se les conoce como argumentos y todas las funciones en R cuentan con al menos un argumento, esto permite personalizar la experiencia del usuario con cada función.

Ahora, si se desea crear un vector repetición bastará con usar la siguiente instrucción.

r1 = rep(1,10) ## Repetir el 1,  10 veces
print(r1)

 [1] 1 1 1 1 1 1 1 1 1 1

r2 = rep(NA,15) ## Repetir NA 15 veces.  
print(r2)

 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

Nota: Los NA se entenderán como valores perdidos. Piense en una base de datos real que presenten vacíos en sus celdas, estos se presentarán como valores NA en R.

Ahora, Observe los siguiente ejemplos

r3 = rep(a, 3)
print(a);print(r3)

[1]  2  4  6  8 10

 [1]  2  4  6  8 10  2  4  6  8 10  2  4  6  8 10

r4 = rep(s, 3)
print(s);print(r4)

[1] "A" "B" "C" "A" "C" "D" "A"

 [1] "A" "B" "C" "A" "C" "D" "A" "A" "B" "C" "A" "C" "D" "A" "A" "B" "C"
[18] "A" "C" "D" "A"

** ¿qué sucede al aplicar la función rep() a un vector?**

Función Ayuda

Como vimos, cada función cuenta con un conjunto de argumentos que serán muy útiles en programación avanzada. Para acceder a los distintos argumentos de la función se puede usar la tecla tab o se usa la función help para una documentacion detallada de la función. Consultemos

help(rep)
?rep  ## equivalente a la funcion help()

** ¿Qué resultado se despliega de la función consultada?**

Nota: Esta función será extremadamente útil cuando se tenga duda de uso de cada función, el entorno de R brindará una documentación detallada de su uso, incluyendo descripción de argumentos y ejemplos.

Instalación de paquetes

Un paquete en R es un conjunto estructurado de funciones y datos que permiten mejorar y complementar los análisis en R. Los paquetes en R pueden ser instalados mediante la instrucción install.packages(). Para instalar un paquete, el nombre deberá ingresarse como argumento de la función entre comillas, como se muestra a continuación

install.packages("ISLR")

Este comando descarga el paquete ISLR y lo instala en el computador. Para hacer uso de las funciones del paquete debera ser cargado en la interfaz y esto se realiza mediante el comando library(), en este caso, las comillas no son necesarias. Para el ejemplo anterior se tiene,

library(ISLR)

Carga de conjuntos de datos desde una librería

Muchos de las librerías que se instalan y cargan en R, contienen conjuntos de datos estructurados (Dataframe) que facilitan el uso de algunas de las funciones y ayudan al usuario a realizar fácilmente cálculos o simulaciones. La librería datasets contiene una gran cantidad de datos con las cuales el usuario puede interactuar, para esto será necesario instalar la librería y usar la instruccionesdata() para obtener el listado de paquetes disponibles.

# install.packages("datasets")
library(datasets)
data()

Para llamar un conjunto de datos específico, llame dentro de la función data() el nombre de los datos entre paréntesis, por ejemplo

library(ISLR)
data("Wage")

Así, se llama el conjunto de datos Wage de la librería ISLR. Estos datos quedan almacenados con el mismo nombre y se puede trabajar con ellos facilmente

names(Wage)  

 [1] "year"       "age"        "maritl"     "race"       "education" 
 [6] "region"     "jobclass"   "health"     "health_ins" "logwage"   
[11] "wage"      

head(Wage)  

       year age           maritl     race       education
231655 2006  18 1. Never Married 1. White    1. < HS Grad
86582  2004  24 1. Never Married 1. White 4. College Grad
161300 2003  45       2. Married 1. White 3. Some College
155159 2003  43       2. Married 3. Asian 4. College Grad
11443  2005  50      4. Divorced 1. White      2. HS Grad
376662 2008  54       2. Married 1. White 4. College Grad
                   region       jobclass         health health_ins
231655 2. Middle Atlantic  1. Industrial      1. <=Good      2. No
86582  2. Middle Atlantic 2. Information 2. >=Very Good      2. No
161300 2. Middle Atlantic  1. Industrial      1. <=Good     1. Yes
155159 2. Middle Atlantic 2. Information 2. >=Very Good     1. Yes
11443  2. Middle Atlantic 2. Information      1. <=Good     1. Yes
376662 2. Middle Atlantic 2. Information 2. >=Very Good     1. Yes
        logwage      wage
231655 4.318063  75.04315
86582  4.255273  70.47602
161300 4.875061 130.98218
155159 5.041393 154.68529
11443  4.318063  75.04315
376662 4.845098 127.11574

** ¿Cual es el resultado de las funciones name() y head() ?**

Dataframe

El trabajo con dataframes en R es aconsejable debido a que en la práctica contamos con varias variables y lo mejor es trabajar con ellos de una manera estructurada. Así como vimos que los conjuntos de datos pueden ser cargados desde una librería, el usuario también podrá crear sus propios conjuntos de datos. Suponga las siguientes variables

X1 = c("A","B","C","D","E","F","G","H")
X2 = c(12,13,15,17,8,12,19,12)
X3 = c("azul","amarillo","azul","verde","amarillo","azul","verde","azul")
X4 = c("F","M","F","F","F","M","M","F")
X5 = c(4.3,4,3.5,2.3,4.5,5,3.2,3.7)
X6 = c(1.56,1.76,1.76,1.32,1.45,1.64,1.78,1.82) 

Datos = data.frame(Nombre = X1, Edad = X2, Equipo = X3, Genero = X4, Nota = X5, Estatura = X6)
Datos

  Nombre Edad   Equipo Genero Nota Estatura
1      A   12     azul      F  4.3     1.56
2      B   13 amarillo      M  4.0     1.76
3      C   15     azul      F  3.5     1.76
4      D   17    verde      F  2.3     1.32
5      E    8 amarillo      F  4.5     1.45
6      F   12     azul      M  5.0     1.64
7      G   19    verde      M  3.2     1.78
8      H   12     azul      F  3.7     1.82

Datos$Equipo  ##acceder a una columna del dataframe con el operador $ 

[1] azul     amarillo azul     verde    amarillo azul     verde    azul    
Levels: amarillo azul verde

Se almacenan los datos en una variable llamada Datos, con esta se podrá realizar distintos análisis, desde el cálculo de estadísticas, consultas, subconjuntos, entre otros.

Consultas básicas en un dataframe

Algunas consultas pueden ser realizadas dentro de un dataframe,

Datos[Datos$Edad < 13,]  # Individuos que tengan menos de 13 años

  Nombre Edad   Equipo Genero Nota Estatura
1      A   12     azul      F  4.3     1.56
5      E    8 amarillo      F  4.5     1.45
6      F   12     azul      M  5.0     1.64
8      H   12     azul      F  3.7     1.82

Datos[Datos$Equipo == "verde",] # Individuos que pertenzzecan al equipo verde

  Nombre Edad Equipo Genero Nota Estatura
4      D   17  verde      F  2.3     1.32
7      G   19  verde      M  3.2     1.78

table(Datos$Genero) ## Numero de hombres y mujeres. table() calcula las frecuencias de las clases F y M

F M 
5 3 

Trabajo con filas y columnas del dataframe

Con [,] accedemos a las posiciones de fila y columna en el conjunto de datos. Por ejemplo, quiero si quiero acceder a la

Datos[1,]  # Fila 1

  Nombre Edad Equipo Genero Nota Estatura
1      A   12   azul      F  4.3     1.56

Datos[,2]  # Columna 2

[1] 12 13 15 17  8 12 19 12

Datos[1,2]  # Fila 1 y Columna 2

[1] 12

Estadísticas

La función summary() nos brindará un resumen adecuado de las variables de un conjunto de datos, considere el ejemplo con el dataframe construido previamente

summary(Datos)

     Nombre       Edad           Equipo  Genero      Nota      
 A      :1   Min.   : 8.0   amarillo:2   F:5    Min.   :2.300  
 B      :1   1st Qu.:12.0   azul    :4   M:3    1st Qu.:3.425  
 C      :1   Median :12.5   verde   :2          Median :3.850  
 D      :1   Mean   :13.5                       Mean   :3.812  
 E      :1   3rd Qu.:15.5                       3rd Qu.:4.350  
 F      :1   Max.   :19.0                       Max.   :5.000  
 (Other):2                                                     
    Estatura    
 Min.   :1.320  
 1st Qu.:1.532  
 Median :1.700  
 Mean   :1.636  
 3rd Qu.:1.765  
 Max.   :1.820  
                

Note que nos calcula medidas de centralidad, dispersión y localización. Pero estas pueden calcularse de manera independiente como sigue, considere el cálculo de los cuantiles

quantile(Datos$Estatura)  ## Cuantiles de la variable predeterminados 

    0%    25%    50%    75%   100% 
1.3200 1.5325 1.7000 1.7650 1.8200 

quantile(Datos$Estatura, probs = c(0.1, 0.6, 0.9))

  10%   60%   90% 
1.411 1.760 1.792 

quantile(Datos$Estatura, probs = c(0.15, 0.3, 0.45))

   15%    30%    45% 
1.4555 1.5680 1.6580 

¿Cuál es la diferencia entre las dos líneas de código ejecutadas?

También podremos calcular

mean(Datos$Estatura)  # Media o promedio

[1] 1.63625

var(Datos$Estatura)  ## Varianza

[1] 0.03222679

sd(Datos$Estatura)  ## Desviaciòn estándar

[1] 0.1795182

sort(Datos$Estatura)  ## Ordena el vector, se puede determinar la mediana

[1] 1.32 1.45 1.56 1.64 1.76 1.76 1.78 1.82

min(Datos$Estatura)   # Mínimo

[1] 1.32

max(Datos$Estatura)   # Máximo

[1] 1.82

IQR(Datos$Estatura)   # Rango Intercuartilico

[1] 0.2325

table(Datos$Estatura)  ## Cuenta la frecuencia de los datos, asì podemos determinar la moda

1.32 1.45 1.56 1.64 1.76 1.78 1.82 
   1    1    1    1    2    1    1 

También se pueden crear funciones, como por ejemplo, la que se crea para determinar la mediana de un vector de datos.

mediana = function(vector){
  l = length(vector)
  vector = sort(vector)
  if(l%%2 == 1) {
    mediana  = vector[round(l/2,0)]
  }else{
    mediana = mean(c(vector[l/2],vector[(l/2)+1]))
  }
  return(paste("La mediana es:",mediana) )
}
mediana(Datos$Estatura)

[1] "La mediana es: 1.7"

prueba = c(1,2,3,4,5,6,7,8)
mediana(prueba)

[1] "La mediana es: 4.5"

Gráficos

Gráficos univariados

Diagrama de Tallos y hojas

stem(Datos$Estatura, scale=2)

  The decimal point is 1 digit(s) to the left of the |

  13 | 2
  14 | 5
  15 | 6
  16 | 4
  17 | 668
  18 | 2

Histograma

data("women") ## Usamos datos del sistema
hist(women$weight, nclass=8, main="Histograma de estaturas", ylab="Frecuencia", xlab="Clases", sub="Fuente: Elaboración propia.", col="Darkgray")

Diagrama de barras

barplot(Datos$Estatura, main="Diagrama de barras - Estatura", xlab="Individuo", ylab="Estatura", sub="Elaboración propia", col = "gray")

Diagrama de torta

data("iris")
conteo = table(iris$Species)
pie(conteo, main="Especies de Iris", sub = "Fuente: RStudio - datos Iris", col=rainbow(3))

Boxplot

boxplot(iris$Sepal.Width, main = "Boxplot - Ancho Sépalo", ylab="Ancho Sépalo", col = "lightgrey", sub="Fuente: Elaboración propia. Datos: RStudio - Iris") 
lines( c(0.8,1.2)  ,  rep(mean(iris$Sepal.Width), 2) , lwd = 2, col=2 )  ## Obligatorio
points(rep(1,length(iris$Sepal.Width)),iris$Sepal.Width, col="darkgreen", cex=0.4, pch =3)   ## Opcional
legend("topright", c("Mediana","Media"), lwd=2 , col=1:2) # Obligatorio

Gráficos bivariados y multivariados

Diagrama Dispersión

data("iris")
plot(iris$Sepal.Length, iris$Sepal.Width, main="Gráfico de dispersión", xlab="Largo Sépalo", ylab="ancho Sépalo", sub="Fuente: Elaboraciòn propia")

Se pueden cambiar algunos elementos de gràfico con los argumentos pch (valores de 0.25), cex y type (“l”, “b”)

plot(iris$Sepal.Length, iris$Sepal.Width, main="Gráfico de dispersión", xlab="Largo Sépalo", ylab="Ancho Sépalo", sub="Fuente: Elaboraciòn propia", pch = 17, cex= 0.8)

plot(1:length(iris$Sepal.Width), iris$Sepal.Width, main="Gráfico de dispersión", xlab="Índice", ylab="ancho Sépalo", sub="Fuente: Elaboraciòn propia", type = "b")

Paquete ggplot2

El paquete ggplot2 (https://ggplot2.tidyverse.org/index.html) nos permitirá generar gráficos sofisticados y es considerado por muchos el paquete de gráficos estadísticos más importantes de R, pertenece al tidyverse que es un conjunto de paquetes especializados para la ciencia de datos (Data Science). El ‘gg’ viene del inglés grammar graphics o gramática de gráficos. A partir de un conjunto de datos, un sistema de coordenadas y un geometrías podremos construir cualquier gráfico. También cuenta con un conjunto de temas que permiten personalizar mucho más nuestros gráficos. Para ver algunos temas y ejemplos ingrese en el siguiente enlace http://www.ggplot2-exts.org/ggthemes.html

Para hacer nuestros primeros gráficos sofisticados, en primer lugar instalamos el paquete ggplot2 mediante el comando install.packages() y lo llamamos mediante el comando library(ggplot2).

install.packages("ggplot2")

library(ggplot2)

Warning: package 'ggplot2' was built under R version 3.5.2

Nota: Cuando el paquete se ha instalado la instrucción install.packages() puede borrarse del script, esto para que no se instale el paquete cada vez que el codigo se ejecuta.

Una vez cargada la librería procedemos con los gráficos básicos sobre el programa. Recuerde que hay elementos básicos cuando construimos un gráfico,

1. la función ggplot() que contiene las variables a gráficar y la función aes()
2. La función de geometría que comienza con geom_*
3. La función labs()
4. (Opcional) funciones de coordenadas, tal como coord_flip()

Graficos de dispersión trivariado

A diferencia del histograma bivariado obtenido mediante la función plot() de la librería base de R, con ggplot2 se podrá realizar un histograma hasta de 5 variables, considerando una variable diferente para el eje x, eje y, el color de acuerdo a una categoría, el tamaño del marcador y la forma. El siguiente es un ejemplo de histograma trivariado con el conjunto de datos iris.

data("iris")
library(ggthemes)
ggplot(iris, aes(x= Sepal.Width, y = Petal.Width, col = Species)) + geom_point() + theme_bw() 

Nota: La librería ggthemes tiene temas adicionales para personalizar gráficos, para el anterior gráfico se usó el tema theme_wsj(), así que si desea usar este estilo u otro debe instalar esta librería con el comando install.packages(ggthemes) y llamarla posteriormente con el comando library(ggthemes). Para ver más información sobre los diferentes temas consulte este enlace https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

data("iris")
library(ggthemes)
ggplot(iris, aes(x= Sepal.Width, y = Petal.Width, shape = Species)) + geom_point() + theme_classic()

data("iris")
library(ggthemes)
ggplot(iris, aes(x= Sepal.Width, y = Petal.Width, size = Species)) + geom_point() + theme_minimal()

Warning: Using size for a discrete variable is not advised.

¿Cuáles son las diferencias, aparte de temas diferentes, de los tres anteriores gráficos?

library(ggthemes)
petalL = cut(iris$Petal.Length, breaks = 3)
ggplot(iris, aes(x= Sepal.Width, y = Petal.Width, col = Species, shape = petalL, size = Sepal.Length)) + geom_point(alpha = 1/3) + theme_bw()

El anterior gráfico es un ejemplo de un gráfico de 5 variables, aunque visualmente no es muy claro es una muestra del poder de la función ggplot(). ¿Puede proponer un mejor gráfico multivariado con iris u otro conjunto de datos?

Gráficos de caja bivariados

Volvamos a los gráficos bivariados que son materia de interés en el curso, para esto supongamos que queremos comparar el ancho del sépalo por especie en el conjunto de datos de iris. Observe el siguiente ejemplo

ggplot(iris, aes(x = Species, y = Sepal.Width)) + geom_boxplot() + labs(x = "Especies" , y = "Ancho Sépalo", title = "Comparación entre especies - Ancho sépalo", subtitle = "Basado en conjunto de datos Iris", caption = "Fuente: Elaboración propia" )

¿Que puede concluir sobre el ancho de sépalo en las tres especies comparadas? ¿qué función cumple la función labs()?

Para consultar información sobre cada geometría use ? por ejemplo,

?geom_point()
?geom_boxplot()

¿qué información obtiene a partir de la ejecución de los anteriores comandos?

Gráficos de barras y barras apiladas

Para ver como realizar gráficos apilados con ggplot2 usaremos el conjunto de datos diamonds, para obtener información sobre este dataset, ejecute el siguiente comando,

?diamonds

¿Qué puede decir de las variables de diamonds?

En los siguientes gráficos se presentan dos ejemplos con las mismas variables, el primer gráfico es un gráfico de barras y el segundo es apilado.

data("diamonds")
ggplot(diamonds, aes(x = cut, fill = color) ) + geom_bar()  + labs(title = "Número de diamantes por corte", subtitle = "Agrupación por color (D-mejor, J-peor)", x = "Corte", y = "Frecuencia", color = "Color", caption= "Fuente: Elaboración propia")

ggplot(diamonds, aes(x = cut, fill = color) ) + geom_bar(position = "fill") + labs(title = "Número de diamantes por corte", subtitle = "Agrupación por color (D-mejor, J-peor)", x = "Corte", y = "Porcentaje", color = "Color", caption= "Fuente: Elaboración propia")

¿Cuál es la diferencia entre las dos líneas de código? Que pasa si queremos hacer la presentación del gráfico con las barras horizontales, agregue la función coord_flip() a alguno de los anteriores gráficos.

En el siguiente ejemplo ya se tiene un gráfico más personalizado, usando algunas de los comandos que hemos aprendido y el tema theme_theeconomist()

data("diamonds")
ggplot(diamonds, aes(x = cut, fill = color) ) + geom_bar()  + labs(title = "Número de diamantes por corte", subtitle = "Agrupación por color (D-mejor, J-peor)", x = "Corte", y = "Frecuencia", color = "Color", caption= "Fuente: Elaboración propia") + coord_flip() + theme_economist() + theme(legend.position = "right")

Si queremos cambiar la leyenda de posición agregue theme(legend.position = "top") o la posición de su preferencia. ¿Cómo luce el anterior gráfico?

Gráficos de área

El último de los gráficos bivariados que trataremos es el gráfico de área, para presentar el ejemplo de este tipo de gráficos continuaremos trabajando con el conjunto de datos diamonds,

ggplot(diamonds, aes(x = price, fill = cut)) + geom_area(stat = "bin", alpha = 0.6) + labs(title = "Precio de diamantes", subtitle = "(Por tipo de corte)", x = "Precio", y = "Frecuencia", fill = "Corte",  caption = "Fuente: Elaboración propia") + theme_bw()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Si quiere más información sobre la función geom_area no olvide consultar la ayuda

?geom_area()

Faceado

Será útil cuando tengamos como objetivo generar un gráfico por cada nivel de una variaable categórica, por ejemplo, tomando como base nuvamente los datos Iris y los tres niveles (que corresponden con las especies), considere el sigueinte ejemplo

ggplot(iris, aes(x= Sepal.Width, y = Sepal.Length)) + geom_point(size = 3, alpha = 1/5) + facet_wrap(~Species) + labs(x = "Ancho Sépalo", y = "Longitud Sépalo", title = "Gráfico de dispersión por especie", caption = "Fuente: Elaboración propia") + geom_smooth()

`geom_smooth()` using method = 'loess' and formula 'y ~ x'

¿Cuál es la principal característica del anterior gráfico?. Consulte más información sobre la función geom_smooth(). Cambie en el anterior gráfico la función geom_smooth() por geom_smooth(method = "lm") ¿Puede explicar este cambio?

?geom_smooth()   

Distribuciones de probabilidad

Las distribuciones de probabilidad nos permite describir como los valores de una variable aleatoria se encuentran distribuidos. Existen distribuciones de probabilidad discretas y continuas, algunas de la más conocidas se presentarán en R, donde se podrán generar números aleatorios, esbozar la función de densidad de la distribución, calcular percentiles y probabilidades. Estas actividades dependerán del prefijo usado para cada distribución: r, d, p y q, por ejemplo pnorm() calcula probabilidades de la distribución normal, qbinom() los percentiles de la distribución binomial, dhyper() valores de la función de densidad de la distribución hipergeométrica, y rexp() generación de números aleatorios de la distribución exponencial.

Distribuciones discretas: Binomial, Hipergeométrica y Poisson

Números aleatorios

Los números aleatorios suelen utilizarse para realizar simulaciones estadística. Para el caso de la distribución binomial se generarán con el comando rbinom() y requerirá del número de intentos (size =) del experimento, la probabilidad de exito (prob =) y el número de realizaciones que se desean generar, considere los siguientes ejemplos,

rbinom(n =  1, prob = 0.2, size = 6)  ### Generar 1 numero aletorio dist binomial con prob 0.2.

[1] 2

rbinom(n =  4, prob = 0.5, size = 6)  ### Generar 4 numero aletorio dist binomial con prob 0.5.

[1] 2 4 5 1

Para el caso de la distribución hipergeométrica, solo bastará con usar el comando rhyper() y la distribución poisson rpois() y siempre se debe definir el tamaño de la muestra a generar y cambiar los parámetros de la distribución. Como son valores aleatorios, estos se generan a partir de los parámetros de las distribucion, es decir, que el conjunto de datos extraidos es una muestra.

Cálculo de probabilidades y Cuantiles

Para el cálculo de probabilidades, solo basta con cambiar la letra inicial del comando por p y para cuantiles por q, en este caso los argumentos de la función para el cálculo de probabilidades son el cuantil en el que se desea calcular la probabilidad y los parámetros de la distribución, por ejemplo para el caso binomial

pbinom(q = 2, prob = 0.2, size = 10, lower.tail = FALSE)

## [1] 0.3222005

pbinom(q = 2, prob = 0.2, size = 10, lower.tail = TRUE)

## [1] 0.6777995

¿Cuál es la diferencia en resultados del lower.tail falso o verdadero ?

Para la distribución hipergeométrica y Poisson se calculan con la función phyper(), ppois(). Para el cálculo de cuantiles la función espera el valor de probabilidad en el cuantil a determinar y los parámetros de la distribución, por ejemplo para el caso binomial.

qbinom(p = 0.322, size = 10, prob = 0.2, lower.tail = FALSE)

## [1] 3

qbinom(p = 0.322, size = 10, prob = 0.2, lower.tail = TRUE)

## [1] 1

¿Cuál es la diferencia en resultados del lower.tail falso o verdadero ?

Función de probabilidad

Usualmente se usa para obtener gráficos de la distribución, o para obtener valores de la densidad en un punto específico y solo basta con cambiar la letra inicial por d al igual que los casos anteriores. Para hacer gráficos de la distribución considere el siguiente ejemplo.

x = seq(0,10)
y = dbinom(x, size = 10, prob = 0.2 )
plot(x,y, type = "h", lwd = 6, main = "Función de probabilidad Bin(10,0.2)", ylab = "p(x)", xlab = "x")

Función de distribución

Se considera que la distribución no es continua, por lo que en este caso se considera la instrucción stepfun() para graficar los saltos,

plot(stepfun(1:10,pbinom(1:11, size = 10, prob = 0.2 )),main="Funcion de distribucion B(10,0.2)")

Distribuciones continuas: Normal y exponencial

La generación de número aleatorios, calculo de cuantiles y probabilidades se hace del mismo modo en distribuciones continuas. Sin embargo, se debe tener en cuenta la continuidad al graficar la función de distribución, a diferencia del caso discreto. Considere el caso de una distribución normal con media \(mean = 10\) y desviación \(sd = 2\). Su función de distribución y de densidad,

Función de distribución Normal

theta = 0:20
plot(theta, pnorm(theta, mean = 10, sd= 2), type = "l", main="Funcion de distribucion N(10,2)")

Función de densidad distribución Normal

theta = seq(0,20, 0.01)
plot(theta, dnorm(theta, mean = 10, sd= 2), type = "l", main="Funcion de densidad N(10,2)")

Usualmente la distribución se puede emplear para visualizar la regiones donde se acumula cierta probabilidad y esto se puede mediante la función polygon().

theta = seq(0,20, 0.01)
plot(theta, dnorm(theta, mean = 10, sd= 2), type = "l", main="Funcion de densidad N(10,2)")

region=seq(10,13,0.01)
xP <-c(10,region,13)
yP <-c(0,dnorm(region,mean = 10,sd = 2),0)
polygon(xP, yP, col = "gray")

Intervalos de confianza y pruebas de hipótesis

Intervalos de confianza para la media

Aunque la estimación puntual en Estadística es muy útil, también los es la estiamción por intervalos dado que nos da una idea del valor que podría tomar el parámetro verdadero con un nivel de confianza y de error determinados. Es necesario recordar que existen dos tipos de errores que son: error Tipo I y error Tipo II y usualmente se controla el más grave de ellos (tipo I). El error tipo I se asocia al nivel de significancia \(\alpha\) y, en la práctica, se suele asumir valores de \(0.5, 0.1, 0.01\) que están muy relacionados con el proceso de muestreo y el campo en el que se esté trabajando, por ejemplo en medicina e industria los valores difieren significativamente en algunas pruebas llegando a sumir valores de \(\alpha = 0.001\). El nivel de confianza de un intervalo se define como \((1-\alpha)\%\) y construimos un intervalo de confianza para la media a partir de la expresión \[ \overline{x} - Z_{\alpha/2} \frac{S }{\sqrt{n}} \leq \mu \leq \overline{x} + Z_{1-\alpha/2} \frac{S}{\sqrt{n}} \]

El comportamiento esperado es que el parámetro verdadero caiga en este intervalo un \((1-\alpha)\%\) de las veces, y mediante simulación esto se puede determinar. En la siguiente imagen, extraemos 100 muestras de una distribución normal con media cero, desviación estándar 1 y obtenemos los intervalos de confianza

¿Cuántos intervalos no contienen la media igual a cero? ¿está relacionado esto con el nivel de significancia?

Ahora, la determinación de los intervalos de confianza en R es sencilla y esto se hace considerando cuantiles de la distribución t-student puesto que cuando n es suficientemente grande su aproximación a la normal es clara, para esto replique la siguiente gráfica variando el tamaño de muestra n en una distribución t student.

¿Para que valor de n la distribución T-student? se paroxima a la distribución normal

Para la construcción de los intervalos de confianza de la media, utilicemos la función CI() del paquete Rmisc y los datos women de la librería datasets.

library(Rmisc) # cargamos la librería Rmisc

## Loading required package: lattice

## Loading required package: plyr

data(women)   # Cargamos los datos 
estatura = women$height
peso = women$weight

## intervalo de confianza para la media 
CI(estatura, ci =  0.95)

##    upper     mean    lower 
## 67.47659 65.00000 62.52341

Funciones en R

Para validar los resultados de la función CI() del paquete Rmisc programemos una función en R. Esto se puede realizar considerando el comando function(), definiendo entre {} los argumentos que trabajaremos dentro de la función y teniendo claro que queremos que la función retorne: número, caracter, dataframe, lista, etc. Para esto, creamos la función que calcula la media de un conjunto de datos, como sigue

mediaDatos = function(datos){
  return(mean(datos))
}
mediaDatos(estatura)

## [1] 65

mean(estatura)

## [1] 65

Usando condicionales if

Ahora creemos la función que nos calcula el cuantil de la distirbución Z o T-Student de acuerdo al tamaño de muestra, supongamos que después de \(n=30\) la distribución T-student se aproxima en distribución a la Normal, para esto usemos el if{}else{} para condicionar sobre el tamaño de muestra, así:

cuantilZT = function(prob, dato){
  if(length(dato)<30){
     cuantil = qt(p = prob, df = (length(dato)-1))
  }else{
     cuantil = qnorm(p = prob)
  }
  return(cuantil)
}
cuantilZT(dato = estatura, prob = 0.95)

## [1] 1.76131

cuantilZT(rep(estatura, 3), prob = 0.95)

## [1] 1.644854

¿Por qué hay diferencias entre el cuantil calculado con la función cuantilZT que creamos?

Ejercicio: Programe una función que calcule el intervalo de confianza para un conjunto de datos con una probabilidad prob y devuelva un arreglo de tres posiciones, limite superior, inferior y la media muestral. Además use el if() para crear un intervalo para proporciones si la variable es de categorica de 2 clases

Pruebas de hipótesis

El contraste de hipótesis está muy relacionado con el cálculo de intervalo de confianza puesto que, de fondo, decidimos sobre una hipótesis planteada de antemano, si el parámetro de interés (e.g. media, proporción) se encuentra o no en el intervalo de confianza a un nivel \(\alpha\) determinado. Para juzgar las pruebas de hipótesis intente tener en cuenta los siguientes pasos para una prueba de hipótesis de media poblacional.

• Defina el nivel de significancia \(\alpha\).
• Postule su hipótesis nula y alternativa, considerando que la nula \(H_o: \mu = \mu_o\) contenga el igual. Hay tres casos posibles de contrastar \[1)Ho: \mu=\mu_o, \hspace{2mm} Ha: \mu\neq\mu_o\] \[2)Ho: \mu=\mu_o, \hspace{2mm} Ha: \mu < \mu_o\] \[3)Ho: \mu=\mu_o, \hspace{2mm} Ha: \mu > \mu_o\]
• Establezca un estadístico de decisión basado en el parámetro de interés y su distribución asociada.
• Defina un criterio de decisión (Región de rechazo) ya sea mediante comparación de cuantiles o comparación de probabilidades (p valor)
• Tome la decisión: rechace o no rechace la hipótesis nula.

En R usamos la función t.test() que además reporta un intervalo de confianza. En la función t.test() los argumentos claves a considerar son el tipo de hipótesis alternativa: alternative = c("two.sided", "less", "greater") y el nivel de confianza: conf.level = 0.95. Así, con lo datos de women() queremos probar la hipótesis a un nivel de confianza del \(95\%\) que la media de la estatura de las mujeres es 60 pulgadas y por lo tanto, formulamos nuestra hipótesis, como sigue \[H_o: \mu = 60, \hspace{2mm} H_a: \mu \neq 60\].

estatura = women$height
t.test(estatura, mu = 60,  alternative = "two.sided", conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  estatura
## t = 4.3301, df = 14, p-value = 0.000692
## alternative hypothesis: true mean is not equal to 60
## 95 percent confidence interval:
##  62.52341 67.47659
## sample estimates:
## mean of x 
##        65

Con un nivel de significancia del \(5\%\) se rechaza la hipótesis nula, es decir, se puede concluir que la estatura media de las mujeres no es de 60 pulgadas (¿Qué pasa si fuese una muestra representativa de la población?).

Si verifica el intervalo de confianza que genera la función t.test() se dará cuenta que el valor de 62 pulgadas de la hipótesis no se encuentra incluido en el intervalo.

Análisis de Varianza ANOVA

A lo largo del curso estuvimos interesados en el análisis de conjunto de variables desde un punto de vista descriptivo e inferencial; esto último con pruebas de hipótesis e intervalos de confianza. Ahora, entendidos los conceptos previos se desea explorar una de las herramientas comunmente usadas en el análisis de la variabilidad de un conjunto determinado de datos para evaluar si existe o no diferencias entre grupos de individuos y conocido como análisis de varianza - ANOVA, donde se analiza la relación de una variable categórica X (con grupos independientes) y sus efectos sobre una variable respuesta Y(numérica). Algunos ejemplos abordados con ANOVA se presentan a continuación:

• Se desea medir la pérdida de peso en un grupo de personas que se someten a 4 dietas diferentes: Control, dieta 1, dieta 2 dieta 3. La variable respuesta es el peso final de la persona y la variable x es el grupo al cual perteneció la persona con una dieta especifica.

• Para medir la efectividad de un medicamento se tienen dos grupos de personas a los cuales se les aplican dos tratamientos: placebo y medicamento de interés. Y se desea evaluar si existen diferencias entre el grupo de placebo y el del medicamento.

La construcción del análisis de varianza ANOVA requiere el cumplimiento de los siguientes supuestos:

• Observaciones independientes.
• Grupos independientes
• Varianza constante (Homocedasticidad).
• Cada grupo tiene un número importante de observaciones (n>20).
• Errores siguen una distribución normal

Cuando n<20 suelen usarse metodología no paramétricas como Kruskal-Wallis

La prueba de hipótesis del análisis de varianza está definida como: \[ Ho: \mu_1 = \mu_2 = \dots = \mu_n \] \[ Ha: \mu_i \neq \mu_j, \hbox{para algún } i \neq j \]

Considere el siguiente ejemplo sobre pesos en seco de unas plantas sometidas a dos tratamientos A y B:

pesos = c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14, 4.81, 4.17,4.41, 3.59, 5.87, 3.83, 6.03,4.89,4.32,4.69, 6.31, 5.12, 5.54, 5.5, 5.37, 5.29, 4.92, 6.15, 5.80, 5.26)
Tratamiento = c(rep("control",10), rep("Tratamiento A",10), rep("Tratamiento B",10))
summary(aov(pesos~Tratamiento))

##             Df Sum Sq Mean Sq F value Pr(>F)  
## Tratamiento  2  3.766  1.8832   4.846 0.0159 *
## Residuals   27 10.492  0.3886                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se puede ver que con un p-valor de 0.0159, a un nivel del \(5\%\), se rechaza la hipótesis nula de igualdad de medias entre tratamientos, esto quiere decir que si hay una diferencia entre los tratamientos usados en las plantas y sus pesos en seco.

Análisis de regresión

El análisis de regresión es una importante técnica estadística para modelizar la relación existente entre una variable respuesta \(y\) (continua) y un conjunto de regresoras o variables explicativas \(x_1, x_2, \dots, x_n\). Por simplicidad, abordaremos el caso en que solo contamos con una regresora \(x_1\), esto es conocido como un análisis de regresión simple y aquí el método de estimación es el mismo que en el caso múltiple: vía mínimos cuadrados ordinarios - MCO. La idea detrás del MCO es minimizar la distancia de cada observación a una recta que ajusta la nube de puntos, la ventaja de MCO es que resulta en estimaciones insesgadas de los parámetros de la recta regresión, es decir \(E(\hat\alpha ) = \alpha\), \(E(\hat\beta ) = \beta\).

El proceso de ajuste de la recta de regresión no solo aborda la estimación, es muy importante la verificación de los supuestos (normalidad y homocedasticidad) sobre los cuales se construye teóricamente el modelo, revisión que se hace sobre los residuales definidos como \(y_i - \hat y_i\) (diferencia entre el valor real y el estimado) y si existen o no valores influyentes y outliers de regresión.

Estimación de la recta de regresión

El método de estimación de los parámetros ampliamente utilizado es el de mínimos cuadrados ordinarios que se obtienen al minimar la distancia de la observación a la recta estimada \[ \sum_i ^n (y_i - \hat y_i)^2\] donde \[\hat y_i = \hat \alpha + \hat \beta x_1\], con \(x_1\) variable explicativa. Por lo tanto, el problema del MCO se reduce a minimizar

\[ \min \bigg\{ \sum_i ^n (y_i - \hat \alpha - \hat \beta x_1)^2 \bigg\}\] Así, derivando respecto a \(\hat\alpha, \hat\beta\) e igualando a cero obtenemos los estimadores insesgados de los parámetros de la regresión, \[\hat \beta = \frac{\sum{(x_i-\overline{x})(y_i-\overline{y})}}{\sum{(x_i-\overline{x})^2}}; \hspace{10mm} \hat\alpha = \overline{y} - \hat\beta\overline{x} \]

Para hacer estimaciones en R consideremos el dataset cars y la función lm() (linear model) definiendo la relación funcional entre las variables con el operador ~. Consideremos la relación entre la velocidad (regresora) y la distancia recorrida (respuesta) y ajustemos el modelo

data(cars)
fit.model = lm(dist ~ speed , data= cars)
fit.model

## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Coefficients:
## (Intercept)        speed  
##     -17.579        3.932

Las estimaciones son $= -17.57 $ y \(\hat \beta = 3.932\). Y si queremos obtener toda la inferencia sobre los parámetros usamos la función summary(), como sigue

summary(fit.model)

## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

En esta salida se nos presenta la estimación, el valor t para juzgar la hipótesis de los coeficientes individuales y el correspondiente p-valor. Las hipótesis que se juzgan en este caso son:

\[ H_o: \alpha = 0, \hspace{10mm} H_a: \alpha \neq 0 \] \[ H_o: \beta = 0, \hspace{10mm} H_a: \beta \neq 0 \]

El p valor nos permite tomar la decisión. En ambos casos a un nivel de significancia \(\alpha = 0.05\) se rechazan ambas hipótesis nulas por lo que se puede concluir que las variables son significativas de manera independiente al modelo de regresión.

Otro elemento importante a considerar es el \(R^2\) o también conocido como coeficiente de determinación que no es más que una indicador de bondad de ajuste del modelo. Este explica el porcentaje de varianza de la variable repuesta que es explicado por la variable explicativa. Es importante mencionar que el \(R^2 \hspace{2mm} ajustado\) no tiene interpretacion en el análisis de regresión simple, pero si es muy importante el \(R^2 \hspace{2mm} múltiple\). Averigue como interpretar el \(R^2\)

Verificación de supuestos

La verificación de los supuestos de la regresión se hace sobre los residuales (errores) y para este nivel diremos que sobre estos recaen supuestos importantes:

*Los errores siguen una distribución normal con media cero y desviación estándar \(\sigma\) (varianza constante).

Aquí haremos una revisión a los residuales de manera gráfica, pues no entraremos en mediciones rigurosas de estos. Por lo tanto, se hará una verificación a manera descriptiva a partir de gráficos ajustados vs residuales estandarizados, histogramas y qq-plot.

Extraemos los residuales y valores ajustados del modelo

residuales = fit.model$residuals
res.estandarizados = scale(residuales, scale = T, center = T)
ajustados = fit.model$fitted.values

Se estandarizan los residuales pues brindan mayor interpretabilidad. Luego generamos el gráfico ajustados vs residuales.std y el histograma

plot(ajustados, res.estandarizados, main = "Ajustados vs residuales", xlab = "Valores ajustados", ylab = "residuales estandarizados" )

hist(res.estandarizados, main = "Histograma de residuales", ylab = "Frecuencias", xlab = "Residuales Estandarizados")

Vemos que algunos valores se alejan más allá de dos (2) desviaciones estándar por lo que se debe revisar su influencia sobre del modelo, posiblemente estas observaciones sean valores atípicos. En el histograma se perciben valores que hacen que la distribución sea sesgada. Para verificarlo generamos el qqplot

qqnorm(res.estandarizados)
abline(0,1, col =2, lwd =2)

En este gráfico se espera que los residuales que se distribuyen normal se agrupen alrededor de la recta de \(45º\), si los valores se alejan mucho da una idea de que el supuesto de normalidad no se está cumpliendo.

Predicciones

El interés sobre muchos modelos de regresión recae sobre su capacidad de hacer predicciones sobre nuevas observaciones. Suponga que queremos conocer la distancia estimada para velocidades de 30 y 32 millas, para esto creamos un dataframe con las nuevas observaciones donde se conserve el mismo nombre de la variable explicativa.

velocidad = data.frame(speed = c(30,32))

y se usa la función predict() para obtener las predicciones sobre el modelo fit.model

predict(fit.model, newdata = velocidad)

##        1        2 
## 100.3932 108.2580

Aquí es importante el cálculo de intervalos de confianza sobre las predicciones, para esto use el argumento interval = "prediction" y el nivel de confianza level = 0.95.

predict(fit.model, newdata = velocidad, interval = "prediction", level = 0.95)

##        fit      lwr      upr
## 1 100.3932 66.86529 133.9210
## 2 108.2580 74.08678 142.4292

También se pueden obtener los intervalos de confianza de la recta de regresión con el argumento interval = "confidence de la misma función y con esto construir las bandas de confianza.

predicciones = predict(fit.model, interval = "confidence")
plot(dist ~speed, data =cars, xlab = "Velocidad", ylab = "Distancia", main = "Bandas de confianza de la recta estimada")
abline(fit.model, col = 2, lwd = 2) ## Recta ajustada
a = lines(cars$speed, predicciones[,2], lwd = 2, col = 4)  ## Límite inferior
b = lines(cars$speed, predicciones[,3], lwd = 2, col = 4)  ## Límite superior

Analice las bandas de confianza de la predicción